################################################################
# Replication file for Braumoeller, Marra, Radice and Bradshaw,
# "Flexible Causal Inference for Political Science,"
# Political Analysis (forthcoming)
################################################################

################################################################
# Generate Figure 1
################################################################

library(copula)
library(VineCopula)
#to produce contour plots of copula

pdf(file="Figure-1.pdf", width=12, height=6)
tau1=2/9
tau2=0.3
tau=0.5
amhMvd <- mvdc(copula = amhCopula(param = iTau(amhCopula(), tau2)),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))


claMvd  <- mvdc(copula = archmCopula(family = "clayton", 
               param = iTau(claytonCopula(), tau) ),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))


fgmMvd  <- mvdc(copula = fgmCopula(param = iTau(fgmCopula(), tau1)),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))


fraMvd <- mvdc(copula = archmCopula(family = "frank", param = iTau(frankCopula(), tau)),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0,sd = 1), 
                              list(mean = 0, sd = 1)))

norMvd  <- mvdc(copula = normalCopula(param = iTau(normalCopula(), tau)),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))
 
gumMvd3 <- mvdc(copula = archmCopula(family = "gumbel", param = iTau(gumbelCopula(), tau)),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))

joeMvd3 <- mvdc(copula = archmCopula(family = "joe", param =iTau(joeCopula(), tau) ),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))

tMvd    <- mvdc(copula = tCopula(param = iTau(tCopula(), tau), df = 3),
               margins = c("norm", "norm"), 
               paramMargins = list(list(mean = 0, sd = 1), 
                              list(mean = 0, sd = 1)))


par(mfrow = c(2, 4), mar = c(2, 2, 1, 1), oma = c(1, 1, 0, 0), mgp = c(2, 1, 0))
contour(amhMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "AMH", cex = 1.3)
contour(claMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Clayton", cex = 1.3)
contour(fgmMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "FGM", cex = 1.3)
contour(fraMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Frank", cex = 1.3)
contour(norMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Gaussian", cex = 1.3)
contour(gumMvd3, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Gumbel", cex = 1.3)
contour(joeMvd3, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Joe", cex = 1.3)
contour(tMvd, dMvdc, xlim = c(-3, 3), ylim = c(-3, 3), box01 = F)
text(0,2.8, "Student-t", cex = 1.3)
dev.off()

################################################################
# Data analysis
################################################################

library(SemiParBIVProbit)
library(foreign)

################################################################
# Load von Stein data from reanalysis by Simmons and Hopkins
################################################################

vonStein.data <- read.csv("apsr2-29-04ab.csv")

##############
# MANSKI BOUND 
##############

mb(vonStein.data$article8, vonStein.data$restrict, Model= "B")

################################################################
# Create variables for temporal splines, signatory status
################################################################

vonStein.data$yrfromj[is.na(vonStein.data$yrfromj)] <- (1900+vonStein.data$year[is.na(vonStein.data$yrfromj)]) - vonStein.data$joinyear[is.na(vonStein.data$yrfromj)]
vonStein.data$yrfromj[vonStein.data$yrfromj<0] <- 0

vonStein.data$lastrest.trunc <- vonStein.data$lastrest
vonStein.data$lastrest.trunc[vonStein.data$lastrest.trunc > 10] <- 10

vonStein.data$YearArt8plus19 <-  (1900+vonStein.data$year) - vonStein.data$YearArt8
vonStein.data$pre.sign <- 0
vonStein.data$pre.sign[vonStein.data$YearArt8plus19 < 0] <- 1

################################################################
# Endogenize Article 8 ratification in a recursive probit
# formulation.
################################################################

eq1.ratify <- article8 ~ useimfcr + openness + chgdp + res_gdp + democracy2 + 
                         gattwto + univers + regnorm + flexible + surveil + 
                         gnpcap + volres + year + s(yrfromj)
                         
eq2.noncomply <- restrict ~ article8 + totvol + bopgdp + useimfcr + openness + 
                            chgdp + res_gdp + democracy2 + gattwto + 
                            as.factor(lastrest.trunc)
                            
theta.eq <- ~ s(YearArt8plus19) #, by = as.factor(pre.sign) )                             

################################################################
# Calculate results for Gaussian (8.1) and AMH (8.2) copulae,
# logit marginals; check to see whether one is preferred
################################################################

result8.1  <- SemiParBIVProbit(list(eq1.ratify, eq2.noncomply, theta.eq), data = vonStein.data, BivD="N", margins=c("logit", "logit"))
result8.2  <- SemiParBIVProbit(list(eq1.ratify, eq2.noncomply, theta.eq), data = vonStein.data, BivD="AMH", margins=c("logit", "logit"))
                            
summary(result8.1)

################################################################
# Generate results for Table 2 (for Table 1 and Appendix
# results, see separate simulation replication file)
################################################################
summary(result8.2)

AIC(result8.1, result8.2)
BIC(result8.1, result8.2)

VuongClarke(result8.1, result8.2)

################################################################
# Check shape of spline
################################################################

plot(result8.2, eq=3, pages=1, scale=0, seWithMean = TRUE, jit= TRUE, shade = TRUE)

################################################################
# Use Bayesian posterior simulation to generate estimates
# of ATE
################################################################

AT(result8.2, eq = 2, nm.end = "article8", hd.plot = TRUE, n.sim = 1000)


names.var <- c("article8", "useimfcr", "openness", "chgdp", 
               "res_gdp", "democracy2", "gattwto", "univers", 
               "regnorm", "flexible", "surveil", "gnpcap", 
               "volres", "year", "restrict", "totvol", "bopgdp",
               "pre.sign", "lastrest.trunc", "YearArt8plus19")  
 
 
# new dataset:
# we do not have missingness on covariates and there is not a mismatch 
# when calculating the ATE

vonStein.data1 <- na.omit(vonStein.data[,names.var]) 
set.seed(1)
 
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==-3 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==-2 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==-1 )
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==0 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==1 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==2 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==3 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==4 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==5 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==6 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==7 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==8 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==9 )  
AT(result8.2, eq = 2, nm.end = "article8", ind = vonStein.data1$YearArt8plus19==10 )  
 
CItau <- round(summary(result8.2)$CItau,2)
CItau[vonStein.data1$YearArt8plus19==-10,]
CItau[vonStein.data1$YearArt8plus19==10,]

AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==-3 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==-2 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==-1 )
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==0 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==1 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==2 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==3 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==4 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==5 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==6 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==7 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==8 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==9 )  
AT(result8.2, eq = 2, nm.end = "article8", type = "univariate", ind = vonStein.data1$YearArt8plus19==10 ) 

################################################################
# Generate Figure 2
################################################################

pdf(file="Figure-2.pdf", width=8, height=6)
par(mfrow=c(1,2))
boxplot(result8.2$tau[vonStein.data1$pre.sign==1,], ylim=c(-0.5,0.5), main=expression(hat(tau)), xlab="Non-signatories")
boxplot(result8.2$tau[vonStein.data1$pre.sign==0,], ylim=c(-0.5,0.5), main=expression(hat(tau)), xlab="Signatories")
dev.off()

################################################################
# Test the screening hypothesis
################################################################
set.seed(1)

nosign.sim <- AT(result8.2, eq = 2, nm.end = "article8", n.sim = 10000, ind = vonStein.data1$pre.sign==1)$sim.AT

sign.sim <- AT(result8.2, eq = 2, nm.end = "article8", n.sim = 10000, ind = vonStein.data1$pre.sign==0)$sim.AT

mean(sign.sim)  # -0.05933592
quantile(sign.sim, 0.025) # -0.08607356 
quantile(sign.sim, 0.975) # -0.03444187
mean(nosign.sim) # -0.1048061
mean(sign.sim-nosign.sim) # 0.04547018
ks.test(sign.sim, nosign.sim)  # D = 0.7858, p-value < 2.2e-16

# Explore impact of theta
sign.sim.uni <- AT(result8.2, eq = 2, nm.end = "article8", type="univariate", n.sim = 10000, ind = vonStein.data1$pre.sign==0)$sim.AT
nosign.sim.uni <- AT(result8.2, eq = 2, nm.end = "article8", type="univariate", n.sim = 10000, ind = vonStein.data1$pre.sign==0)$sim.AT
ks.test(sign.sim, sign.sim.uni)  # D = 0.0905, p-value < 2.2e-16
ks.test(nosign.sim, nosign.sim.uni) # # D = 0.8302, p-value < 2.2e-16

mean(sign.sim.uni)  # -0.05688933
quantile(sign.sim.uni, 0.025) # -0.08116282 
quantile(sign.sim.uni, 0.975) # -0.03541766 
mean(nosign.sim.uni)  # -0.05694538
quantile(nosign.sim.uni, 0.025) # -0.08151121 
quantile(nosign.sim.uni, 0.975) # -0.03498611

################################################################
# Generate Figure 3
################################################################

set.seed(1)

full.AT <- AT(result8.2, eq = 2, nm.end = "article8", n.sim=1000, type = "simultaneous", ind = vonStein.data1$pre.sign==1)

full.sign.AT <- AT(result8.2, eq = 2, nm.end = "article8", n.sim = 10000, type = "simultaneous", ind = vonStein.data1$pre.sign==0)


plotmat <- matrix(, nrow = 11, ncol = 3)
plotmat[1,] <- ((AT(result8.2, eq = 2, nm.end = "article8", n.sim=1000,ind = (vonStein.data1$pre.sign==0 & vonStein.data1$YearArt8plus19==0))$res)*100) * -1

for(row in 1:10){
  plotmat[(row + 1) ,] <- ((AT(result8.2, eq = 2, nm.end = "article8", n.sim=1000,ind = vonStein.data1$YearArt8plus19==row)$res)*100) * -1
}

pdf(file="Figure-3.pdf", width=8, height=5)
plot(NULL, xlim=c(-5,9.5), ylim=c(-2,40), ylab=" ", xlab=" ", main="Effects of Article VIII Commitment", xaxt="n", cex.axis=0.7)
axis(1, at=c(0,2,4,6,8,10), cex.axis=0.8)

rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col="#E6E6E6")
abline(h=c(0,10,20,30,40), col="white")
abline(v=c(0,2,4,6,8,10), col="white")

polygon(c(0:10,10:0), c(predict(loess(plotmat[,1]~c(0:10))), predict(loess(rev(plotmat[,3])~c(10:0)))),  border="#00000000", col="#00000033")
lines(c(0:10), predict(loess(plotmat[,2] ~ c(0:10))), col="grey40")
# abline(h=0, lty=2, lwd=0.5)

# polygon(c(0:10,10:0), c(plotmat[,1], rev(plotmat[,3])), border="grey80", col="grey80")
# lines(c(0:10), plotmat[,2])
# abline(h=0, lty=2)

dens.full <- density(full.AT$sim.AT*-100)
polygon(c((dens.full$y*15)-5.6, rep(-5.6, length(dens.full$y))), c(dens.full$x, rev(dens.full$x)), border=NA, col="#00000033")

dens.full.sign <- density(full.sign.AT$sim.AT*-100)
polygon(c((dens.full.sign$y*15)-5.6, rep(-5.6, length(dens.full.sign$y))), c(dens.full.sign$x, rev(dens.full.sign$x)), border=NA, col="#00000033")

# segments(-5, full.AT$res[2]*-100, 0, full.AT$res[2]*-100, lty=3, col="grey50")

text(-3, 36, "Screening", cex=0.7)
text(-3, 34, "effect", cex=0.7)
text(5, 36, "Constraining", cex=0.7)
text(5, 34, "effect", cex=0.7)
mtext("Treatment effect", side=2, line=2, cex=0.7)
mtext("                                                Years since signing", side=1, line=2, cex=0.7)

abline(v=0)
box(lwd=1)

dev.off()



